#---------------------------------#
# Modeling the network diffusion  #
#---------------------------------#
rm(list = ls())

#dir.create("data_visualization", showWarnings = F, recursive = T)
dir.create("model_output", showWarnings = F, recursive = T)
range01 <- function(x){(x-min(x, na.rm = T))/(max(x, na.rm = T)-min(x, na.rm = T))}

## Load in data
adjacency_matricies <- readRDS("intl_order/adjacency_matricies_replication.rds")
polity.diffusion <- readRDS("data/polity_diffusion_replication.rds")
min_year <- min(polity.diffusion$year) - 1


## Create variables
polity.diffusion <- mutate(
  polity.diffusion,
  democ = ifelse(polity >= 6, 1, 0),
  autoc = ifelse(polity <= -6, 1, 0),
  v2csgender_ord = as.numeric(as.character(v2csgender_ord)),
  v2csrlgrep_ord = as.numeric(as.character(v2csrlgrep_ord)),
  parcomp = as.numeric(as.character(parcomp)),
  parreg = as.numeric(as.character(parreg)),
  xconst = as.numeric(as.character(xconst)),
  v2mecenefm_ord = as.numeric(as.character(v2mecenefm_ord)),
  women_suffrage = ifelse(v2csgender_ord >= 3, 1, 0),
  religious_freedom = ifelse(v2csrlgrep_ord >= 3, 1, 0),
  free_media = ifelse(v2mecenefm_ord >= 3, 1, 0),
  wdi_gdpcapgr = as.numeric(scale(wdi_gdpcapgr)),
  comptetive_political_participation = ifelse(parcomp == 5, 1, 0),
  unilateral_executive = ifelse(xconst  == 1, 1, 0),
  restricted_political_participation = ifelse(parreg == 4, 1, 0),
  executive_constrainted = ifelse(xconst == 7, 1, 0)
)


set.seed(1)
training_rows <- sample(nrow(polity.diffusion), 0.9 * nrow(polity.diffusion))
polity.diffusion$year_num <- as.numeric(scale(polity.diffusion$year))
polity.diffusion$gle_cgdpc <- as.numeric(scale(log(polity.diffusion$gle_cgdpc)))
polity.diffusion$wdi_gdpcapgr <- as.numeric(scale(log(polity.diffusion$wdi_gdpcapgr + -1 * min(polity.diffusion$wdi_gdpcapgr, na.rm = T) + 1)))
polity.diffusion.train <- polity.diffusion[training_rows,]
polity.diffusion.test  <- polity.diffusion[-training_rows,]

fe.democ <- fixest::feglm(democ ~ prop_democ_mimic + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war, 
                          data = polity.diffusion.train, family = "binomial")

fe.autoc <- fixest::feglm(autoc ~ prop_autoc_mimic + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war, 
                          data = polity.diffusion.train, family = "binomial")

## Institutional diffusion of political participation restrictions
fe.rest.part <- fixest::feglm(restricted_political_participation ~ prop_rest_pol_part + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war, 
                              data = polity.diffusion.train, family = "binomial")

## Institutional diffusion of free media
fe.free.media <- fixest::feglm(free_media ~ prop_gov_censorship_mimic + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war, 
                               data = polity.diffusion.train, family = "binomial")

## Institutional diffusion of executive constraints
fe.exec.const <- fixest::feglm(executive_constrainted ~ prop_exec_constraint + exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war,
                               data = polity.diffusion.train, family = "binomial")

reduced.democ <- fixest::feglm(democ ~ democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war, 
                          data = polity.diffusion.train, family = "binomial")

reduced.autoc <- fixest::feglm(autoc ~ autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war, 
                          data = polity.diffusion.train, family = "binomial")

## Institutional diffusion of political participation restrictions
reduced.rest.part <- fixest::feglm(restricted_political_participation ~ rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war, 
                              data = polity.diffusion.train, family = "binomial")

## Institutional diffusion of free media
reduced.free.media <- fixest::feglm(free_media ~ gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war, 
                               data = polity.diffusion.train, family = "binomial")

## Institutional diffusion of executive constraints
reduced.exec.const <- fixest::feglm(executive_constrainted ~ exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war, 
                               data = polity.diffusion.train, family = "binomial")




## Institutional diffusion of democracy
suppressWarnings({
  re.democ <- glmer(democ ~ prop_democ_mimic + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                  family = binomial(link = "logit"), 
                  control = glmerControl(
                    optimizer = "bobyqa",
                                       optCtrl = list(maxfun = 5e10), 
                    check.scaleX = "ignore",
                    check.conv.hess = "ignore",
                    check.conv.grad = "ignore"), 
                  data = polity.diffusion.train)
  
  ## Institutional diffusion of autocracy
  re.autoc <- glmer(autoc ~ prop_autoc_mimic + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                  family = binomial(link = "logit"), 
                  control = glmerControl(optimizer = "bobyqa",
                                       optCtrl = list(maxfun = 5e10))
                  , 
                  data = polity.diffusion.train)
  
  
  ## Institutional diffusion of political participation restrictions
  re.rest.part <- glmer(restricted_political_participation ~ prop_rest_pol_part + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                      family = binomial(link = "logit"), 
                      control = glmerControl(optimizer = "bobyqa",
                                             optCtrl = list(maxfun = 5e10))
                      ,  
                      data = polity.diffusion.train)
  
  ## Institutional diffusion of free media
  re.free.media <- glmer(free_media ~ prop_gov_censorship_mimic + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2) + (1|year) + (1|names), 
                       family = binomial(link = "logit"), 
                       control = glmerControl(optimizer = "bobyqa",
                                              optCtrl = list(maxfun = 5e10))
                       ,  
                       data = polity.diffusion.train)
  
  
  ## Institutional diffusion of executive constraints
  re.exec.const <- glmer(executive_constrainted ~ prop_exec_constraint + exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                       family = binomial(link = "logit"), 
                       control = glmerControl(
                         optimizer = "bobyqa", optCtrl = list(maxfun = 5e10)
                       ), 
                       data = polity.diffusion.train)
  
  
})

summary(re.democ)
summary(re.autoc)
summary(re.rest.part)
summary(re.free.media)
summary(re.exec.const)

pred_data <- data.frame(
  prop_exec_constraint = c(0, 1),
  exec_constraint_neighbor = median(polity.diffusion$exec_constraint_neighbor, na.rm = T), 
  
  prop_rest_pol_part = c(0, 1),
  rest_pol_part_neighbor = median(polity.diffusion$rest_pol_part_neighbor, na.rm = T),
  
  prop_gov_censorship_mimic = c(0, 1),
  gov_censorship_neighbor = median(polity.diffusion$gov_censorship_neighbor, na.rm = T), 
  
  prop_autoc_mimic = c(0, 1),
  autocratic_neighbor = median(polity.diffusion$autocratic_neighbor, na.rm = T), 
  
  prop_democ_mimic = c(0, 1),
  democratic_neighbor = median(polity.diffusion$democratic_neighbor, na.rm = T), 
  
  gle_cgdpc = median(polity.diffusion$gle_cgdpc, na.rm = T),
  wdi_gdpcapgr = median(polity.diffusion$wdi_gdpcapgr, na.rm = T),
  civil_war = median(polity.diffusion$civil_war, na.rm = T), 
  year_num  = 0.7181701, 
  names = "MEX",
  year = 1990
)

polity.diffusion.test <- bind_rows(
  polity.diffusion, 
  pred_data
)

democ.pred      <- predict(re.democ, newdata = polity.diffusion.test[(nrow(polity.diffusion.test) - 1):nrow(polity.diffusion.test),], type = "response",  allow.new.levels = T)
autoc.pred      <- predict(re.autoc, newdata = polity.diffusion.test[(nrow(polity.diffusion.test) - 1):nrow(polity.diffusion.test),], type = "response",  allow.new.levels = T)
rest.part.pred  <- predict(re.rest.part, newdata = polity.diffusion.test[(nrow(polity.diffusion.test) - 1):nrow(polity.diffusion.test),], type = "response",  allow.new.levels = T)
free.media.pred <- predict(re.free.media, newdata = polity.diffusion.test[(nrow(polity.diffusion.test) - 1):nrow(polity.diffusion.test),], type = "response",  allow.new.levels = T)
exec.const.pred <- predict(re.exec.const, newdata = polity.diffusion.test[(nrow(polity.diffusion.test) - 1):nrow(polity.diffusion.test),], type = "response",  allow.new.levels = T)

round(exec.const.pred[2] - exec.const.pred[1], 2)
round(rest.part.pred[2] - rest.part.pred[1], 2)
round(free.media.pred[2] - free.media.pred[1], 2)
round(autoc.pred[2] - autoc.pred[1], 2)
round(democ.pred[2] - democ.pred[1], 2)

dir.create("model_output", showWarnings = F, recursive = T)
save(re.democ, file = "model_output/re.democ.rda")
save(re.autoc, file = "model_output/re.autoc.rda")
save(re.rest.part, file = "model_output/re.rest.part.rda")
save(re.free.media, file = "model_output/re.free.media.rda")
save(re.exec.const, file = "model_output/re.exec.const.rda")

## Compare predictive accuracy with and without network ties
full_auc_democ         <- c()
full_auc_autoc         <- c()
full_auc_rest_part     <- c()
full_auc_free_media    <- c()
full_auc_exec_const    <- c()
reduced_auc_democ      <- c()
reduced_auc_autoc      <- c()
reduced_auc_rest_part  <- c()
reduced_auc_free_media <- c()
reduced_auc_exec_const <- c()

set.seed(1)
for (i in 1:5000)
{
  cat("\r")
  cat(paste0("\rPct of boot AUC: ", round(i/5000, 2)))
  boot_test <- polity.diffusion.train[sample(nrow(polity.diffusion.train), nrow(polity.diffusion.train), replace = T),]
  
  full_fit_democ      <- predict(fe.democ, newdata = boot_test, type = "response")
  full_fit_autoc      <- predict(fe.autoc, newdata = boot_test, type = "response")
  full_fit_rest_part  <- predict(fe.rest.part, newdata = boot_test, type = "response")
  full_fit_free_media <- predict(fe.free.media, newdata = boot_test, type = "response")
  full_fit_exec_const <- predict(fe.exec.const, newdata = boot_test, type = "response")
  
  
  reduced_fit_democ      <- predict(reduced.democ, newdata = boot_test, type = "response")
  reduced_fit_autoc      <- predict(reduced.autoc, newdata = boot_test, type = "response")
  reduced_fit_rest_part  <- predict(reduced.rest.part, newdata = boot_test, type = "response")
  reduced_fit_free_media <- predict(reduced.free.media, newdata = boot_test, type = "response")
  reduced_fit_exec_const <- predict(reduced.exec.const, newdata = boot_test, type = "response")
  suppressMessages({
    full_auc_democ[i]         <- as.numeric(pROC::auc(boot_test$democ, full_fit_democ))
    full_auc_autoc[i]         <- as.numeric(pROC::auc(boot_test$autoc, full_fit_autoc))
    full_auc_rest_part[i]     <- as.numeric(pROC::auc(boot_test$restricted_political_participation, full_fit_rest_part))
    full_auc_free_media[i]    <- as.numeric(pROC::auc(boot_test$free_media, full_fit_free_media))
    full_auc_exec_const[i]    <- as.numeric(pROC::auc(boot_test$executive_constrainted, full_fit_exec_const))
    reduced_auc_democ[i]      <- as.numeric(pROC::auc(boot_test$democ, reduced_fit_democ))
    reduced_auc_autoc[i]      <- as.numeric(pROC::auc(boot_test$autoc, reduced_fit_autoc))
    reduced_auc_rest_part[i]  <- as.numeric(pROC::auc(boot_test$restricted_political_participation, reduced_fit_rest_part))
    reduced_auc_free_media[i] <- as.numeric(pROC::auc(boot_test$free_media, reduced_fit_free_media))
    reduced_auc_exec_const[i] <- as.numeric(pROC::auc(boot_test$executive_constrainted, reduced_fit_exec_const))
  })
}


fit_democ      <- predict(re.democ, newdata = polity.diffusion.test, type = "response",  allow.new.levels = T)
fit_autoc      <- predict(re.autoc, newdata = polity.diffusion.test, type = "response",  allow.new.levels = T)
fit_rest_part  <- predict(re.rest.part, newdata = polity.diffusion.test, type = "response",  allow.new.levels = T)
fit_free_media <- predict(re.free.media, newdata = polity.diffusion.test, type = "response",  allow.new.levels = T)
fit_exec_const <- predict(re.exec.const, newdata = polity.diffusion.test, type = "response",  allow.new.levels = T)

table(round(fit_democ), polity.diffusion.test$democ)
table(round(fit_autoc), polity.diffusion.test$autoc)
table(round(fit_rest_part), polity.diffusion.test$restricted_political_participation)
table(round(fit_free_media), polity.diffusion.test$free_media)
table(round(fit_exec_const), polity.diffusion.test$executive_constrainted)

round(sum(diag(prop.table(table(round(fit_democ), polity.diffusion.test$democ)))), 2)
round(sum(diag(prop.table(table(round(fit_autoc), polity.diffusion.test$autoc)))), 2)
round(sum(diag(prop.table(table(round(fit_rest_part), polity.diffusion.test$restricted_political_participation)))), 2)
round(sum(diag(prop.table(table(round(fit_free_media), polity.diffusion.test$free_media)))), 2)
round(sum(diag(prop.table(table(round(fit_exec_const), polity.diffusion.test$executive_constrainted)))), 2)

tp_democ <- table(round(fit_democ), polity.diffusion.test$democ)[2,2]
tp_autoc <- table(round(fit_autoc), polity.diffusion.test$autoc)[2,2]
tp_restp <- table(round(fit_rest_part), polity.diffusion.test$restricted_political_participation)[2,2]
tp_fmedi <- table(round(fit_free_media), polity.diffusion.test$free_media)[2,2]
tp_execc <- table(round(fit_exec_const), polity.diffusion.test$executive_constrainted)[2,2]

tn_democ <- table(round(fit_democ), polity.diffusion.test$democ)[1,1]
tn_autoc <- table(round(fit_autoc), polity.diffusion.test$autoc)[1,1]
tn_restp <- table(round(fit_rest_part), polity.diffusion.test$restricted_political_participation)[1,1]
tn_fmedi <- table(round(fit_free_media), polity.diffusion.test$free_media)[1,1]
tn_execc <- table(round(fit_exec_const), polity.diffusion.test$executive_constrainted)[1,1]

fp_democ <- table(round(fit_democ), polity.diffusion.test$democ)[2,1]
fp_autoc <- table(round(fit_autoc), polity.diffusion.test$autoc)[2,1]
fp_restp <- table(round(fit_rest_part), polity.diffusion.test$restricted_political_participation)[2,1]
fp_fmedi <- table(round(fit_free_media), polity.diffusion.test$free_media)[2,1]
fp_execc <- table(round(fit_exec_const), polity.diffusion.test$executive_constrainted)[2,1]

fn_democ <- table(round(fit_democ), polity.diffusion.test$democ)[1,2]
fn_autoc <- table(round(fit_autoc), polity.diffusion.test$autoc)[1,2]
fn_restp <- table(round(fit_rest_part), polity.diffusion.test$restricted_political_participation)[1,2]
fn_fmedi <- table(round(fit_free_media), polity.diffusion.test$free_media)[1,2]
fn_execc <- table(round(fit_exec_const), polity.diffusion.test$executive_constrainted)[1,2]

precision_democ <- tp_democ / (tp_democ + fp_democ)
precision_autoc <- tp_autoc / (tp_autoc + fp_autoc)
precision_restp <- tp_restp / (tp_restp + fp_restp)
precision_fmedi <- tp_fmedi / (tp_fmedi + fp_fmedi)
precision_execc <- tp_execc / (tp_execc + fp_execc)

recall_democ <- tp_democ / (tp_democ + fn_democ)
recall_autoc <- tp_autoc / (tp_autoc + fn_autoc)
recall_restp <- tp_restp / (tp_restp + fn_restp)
recall_fmedi <- tp_fmedi / (tp_fmedi + fn_fmedi)
recall_execc <- tp_execc / (tp_execc + fn_execc)

f1_democ <- round((2 * precision_democ * recall_democ)/(precision_democ + recall_democ), 2)
f1_autoc <- round((2 * precision_autoc * recall_autoc)/(precision_autoc + recall_autoc), 2)
f1_restp <- round((2 * precision_restp * recall_restp)/(precision_restp + recall_restp), 2)
f1_fmedi <- round((2 * precision_fmedi * recall_fmedi)/(precision_fmedi + recall_fmedi), 2)
f1_execc <- round((2 * precision_execc * recall_execc)/(precision_execc + recall_execc), 2)

boot_auc <- data.frame(
  full_auc_democ,
  reduced_auc_democ,
  full_auc_autoc,
  reduced_auc_autoc,
  full_auc_rest_part,
  reduced_auc_rest_part,
  full_auc_free_media,
  reduced_auc_free_media,
  full_auc_exec_const,
  reduced_auc_exec_const
)

t.test.auc.democ      <- t.test(full_auc_democ, reduced_auc_democ)
t.test.auc.autoc      <- t.test(full_auc_autoc, reduced_auc_autoc)
t.test.auc.rest_part  <- t.test(full_auc_rest_part, reduced_auc_rest_part)
t.test.auc.free_media <- t.test(full_auc_free_media, reduced_auc_free_media)
t.test.auc.exec_const <- t.test(full_auc_exec_const, reduced_auc_exec_const)

saveRDS(boot_auc, file = "model_output/boot_auc_values.rds")
saveRDS(t.test.auc.democ, file = "model_output/t.test.auc.democ.rds")
saveRDS(t.test.auc.autoc, file = "model_output/t.test.auc.autoc.rds")
saveRDS(t.test.auc.rest_part, file = "model_output/t.test.auc.rest_part.rds")
saveRDS(t.test.auc.free_media, file = "model_output/t.test.auc.free_media.rds")
saveRDS(t.test.auc.exec_const, file = "model_output/t.test.auc.exec_const.rds")

pred_democ      <- predict(re.democ, newdata = polity.diffusion.test, type = "response",  allow.new.levels = T)
pred_autoc      <- predict(re.autoc, newdata = polity.diffusion.test, type = "response",  allow.new.levels = T)
pred_rest_part  <- predict(re.rest.part, newdata = polity.diffusion.test, type = "response",  allow.new.levels = T)
pred_free_media <- predict(re.free.media, newdata = polity.diffusion.test, type = "response",  allow.new.levels = T)
pred_exec_const <- predict(re.exec.const, newdata = polity.diffusion.test, type = "response",  allow.new.levels = T)

plot_roc <- data.frame(
  pred = c(
    pred_democ,
    pred_autoc,
    pred_rest_part,
    pred_free_media,
    pred_exec_const
  ),
  result = c(
    polity.diffusion.test$democ,
    polity.diffusion.test$autoc,
    polity.diffusion.test$restricted_political_participation ,
    polity.diffusion.test$free_media,
    polity.diffusion.test$executive_constrainted
  ),
  dv = c(
    rep("Democracy", length(polity.diffusion.test$democ)),
    rep("Autocracy", length(polity.diffusion.test$democ)),
    rep("Restrictions\non political\nparticipation", length(polity.diffusion.test$democ)),
    rep("Free\nmedia", length(polity.diffusion.test$democ)),
    rep("Executive\nconstraints", length(polity.diffusion.test$democ))
  )
) %>% filter(
  !is.na(pred)
) |> mutate(
  facet_lab = "Institutional diffusion"
)

#define object to plot
insti_roc_curve <- ggplot(plot_roc, aes(m = pred, d = result, col = dv)) + 
  geom_roc(n.cuts = 50, labels = F, size = 1) + 
  facet_wrap(~facet_lab) + 
  theme_bw() + 
  scale_colour_brewer(palette = "Dark2") + 
  theme(legend.position = "bottom",
        legend.title = element_blank()) + 
  labs(x = "False positive fraction",
       y = "True positive fraction") + 
theme(
       legend.text = element_text(size = 10),
       strip.text.x = element_text(size = 12, color = "black"), 
       axis.text = element_text(size = 12,  color = "black"), 
       axis.title = element_text(size = 12,  color = "black")) + 
  guides(col = guide_legend(nrow = 2, byrow = TRUE)) +
  geom_abline(slope = 1, intercept = 0, color = "black", size = 1, linetype = "dashed")

t.test.plot <- data.frame(
  lb = c(
    t.test.auc.democ$conf.int[1],
    t.test.auc.democ$conf.int[1],
    t.test.auc.rest_part$conf.int[1],
    t.test.auc.free_media$conf.int[1],
    t.test.auc.exec_const$conf.int[1]   
  ),
  ub = c(
    t.test.auc.democ$conf.int[2],
    t.test.auc.democ$conf.int[2],
    t.test.auc.rest_part$conf.int[2],
    t.test.auc.free_media$conf.int[2],
    t.test.auc.exec_const$conf.int[2]
  ),
  boot = c(
    "Democracy",
    "Autocracy",
    "Restrictions\non political\nparticipation",
    "Free\nmedia",
    "Executive\nconstraints"
  ),
  facet_lab = "Bootstrapped difference full versus restricted AUC"
)


t_test_means <- ggplot(t.test.plot, aes(xmin = lb, xmax = ub, y = reorder(boot, lb), col = boot)) + 
  geom_errorbar(size = 1.5, width = 0.25) + 
  facet_wrap(~facet_lab) + 
  theme_bw() + 
  scale_colour_brewer(palette = "Dark2") + 
  theme(legend.position = "bottom",
        legend.title = element_blank()) + 
  labs(x = "Difference in bootstrapped AUC values (5000 simulations)",
       y = "")+
  theme(legend.position = "none",
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 12, color = "black"), 
        axis.text = element_text(size = 12,  color = "black"), 
        axis.title = element_text(size = 12,  color = "black")) + 
  geom_vline(xintercept = 0, size = 1.5, linetype = "dashed", col = "black")


plot_auc_data <- plot_grid(insti_roc_curve, t_test_means, ncol = 2, align = "h", axis = "bt", rel_widths = c(0.85, 1))
ggsave("data_visualization/diffusion_out_of_sample.png", plot_auc_data, width = 12, height = 7, units = "in")  

set.seed(1)
## Permutation test
polity.diffusion.merge <- dplyr::select(
  polity.diffusion,
  c(
    democ, autoc, executive_constrainted, free_media , restricted_political_participation,
    democratic_neighbor, autocratic_neighbor, exec_constraint_neighbor, gov_censorship_neighbor, rest_pol_part_neighbor,
    gle_cgdpc, wdi_gdpcapgr, civil_war, year_num, year, names,
    v2mecenefm_ord
  )
)
permutation_networks <- list()
X <- 1:500
for (permute_nets in X) 
{
 polity.diffusion.permutation <- list()
 for(j in 1:length(adjacency_matricies)) 
 {
   
   orig_graph    <- reshape2::melt(adjacency_matricies[[j]])
   orig_graph <- subset(orig_graph, Var1 != Var2)
   orig_graph <- mutate(
     orig_graph, 
     Var1  = as.character(Var1),
     Var2  = as.character(Var2),
     Var1a = if_else(Var1 < Var2, Var1, Var2),
     Var2a = if_else(Var2 < Var1, Var1, Var2),
     Var1 = Var1a, 
     Var2 = Var2a
    ) %>% dplyr::select(
      Var1, Var2, value
    ) %>% distinct(
      Var1, Var2, .keep_all = T
    )
   orig_graph$value <- sample(orig_graph$value)
   new_graph <- reshape2::dcast(orig_graph, Var1 ~ Var2)
   rownames(new_graph) <- new_graph$Var1
   new_graph <- new_graph[,-1]
   new_graph <- as.matrix(new_graph)
   g <- graph_from_adjacency_matrix(new_graph, mode = "undirected", weighted = T)
   indeg <- data.frame(
     names = names(degree(g)),
     indegree = strength(g, mode = "in")
   )
   
   
   democ_autoc_data <- filter(
     polity.diffusion, 
     year == j + min_year
   ) %>% 
     distinct(
       names, year, .keep_all = T
     )
   
   power_temp <- distinct(
       democ_autoc_data,
       names, cinc
     )
   
   weighted_edges <- dplyr::rename(
     orig_graph, 
     to = Var1, 
     from = Var2,
     tie = value
   )
   
   edge.indeg <- left_join(
     weighted_edges,
     power_temp,
     by = c("from" = "names")
   ) %>% left_join(
     power_temp,
     by = c("to" = "names")
   ) %>% mutate(
     democ_flag = if_else(from %in% subset(democ_autoc_data, polity >= 6)$names, 1, 0),
     democ_power_flag = if_else(from %in% subset(democ_autoc_data, polity >= 6)$names & cinc.x > cinc.y, 1, 0),
     democ_flag_vdem = if_else(from %in% subset(democ_autoc_data, v2x_polyarchy >= 0.5)$names, 1, 0),
     democ_power_flag_vdem = if_else(from %in% subset(democ_autoc_data, v2x_polyarchy >= 0.5)$names & cinc.x > cinc.y, 1, 0),
     autoc_flag = if_else(from %in% subset(democ_autoc_data, polity <= -6)$names, 1, 0),
     autoc_power_flag = if_else(from %in% subset(democ_autoc_data, polity <= -6)$names & cinc.x > cinc.y, 1, 0),
     exec_constraint = if_else(from %in% subset(democ_autoc_data,  xconst  == 7)$names & cinc.x > cinc.y, 1, 0),
     comp_pol_part = if_else(from %in% subset(democ_autoc_data,  parcomp  == 5)$names & cinc.x > cinc.y, 1, 0),
     unli_exec = if_else(from %in% subset(democ_autoc_data,  xconst  == 1)$names & cinc.x > cinc.y, 1, 0),
     rest_pol_part = if_else(from %in% subset(democ_autoc_data,  parreg  == 4)$names, 1, 0),
     women_participation_flag = if_else(from %in% subset(democ_autoc_data, v2csgender_ord >= 3)$names & cinc.x > cinc.y, 1, 0),
     gov_censorship_flag = if_else(from %in% subset(democ_autoc_data, v2mecenefm_ord >= 3)$names & cinc.x > cinc.y, 1, 0),
     religious_freedom_flag = if_else(from %in% subset(democ_autoc_data, v2csrlgrep_ord >= 3)$names & cinc.x > cinc.y, 1, 0)
   ) %>% plyr::ddply(
     ~to,
     summarize, 
     indeg_democracy = sum(democ_flag * tie),
     indeg_democ_power_int = sum(democ_power_flag  * tie),
     indeg_democracy_vdem = sum(democ_flag_vdem  * tie),
     indeg_democ_power_int_vdem = sum(democ_power_flag_vdem * tie),
     indeg_women_participation = sum(women_participation_flag * tie),
     indeg_autoc_flag = sum(autoc_flag * tie),
     indeg_autoc_power_flag = sum(autoc_power_flag * tie),
     indeg_exec_constraint = sum(exec_constraint * tie),
     indeg_comp_pol_part = sum(comp_pol_part * tie),
     indeg_unli_exec = sum(unli_exec * tie),
     indeg_rest_pol_part = sum(rest_pol_part * tie),
     indeg_gov_censorship = sum(gov_censorship_flag * tie),
     indeg_religious_freedom = sum(religious_freedom_flag * tie)
   ) %>% left_join(
     indeg, 
     by = c(
       "to" = "names"
     )
   )
   edge.indeg <- edge.indeg[!duplicated(edge.indeg),]
   democ_autoc_data <- data.frame(democ_autoc_data)
   democ_autoc_data <- democ_autoc_data[,!(colnames(democ_autoc_data) %in% c(colnames(edge.indeg), colnames(indeg))[c(colnames(edge.indeg), colnames(indeg)) != "names"])]
   suppressWarnings(year.data <- left_join(
     edge.indeg, 
     democ_autoc_data,
     by = c("to" = "names")
   )  %>% mutate(
     indeg_democracy = tidyr::replace_na(indeg_democracy, 0),
     indeg_democ_power_int = tidyr::replace_na(indeg_democ_power_int, 0),
     indeg_democracy_vdem = tidyr::replace_na(indeg_democracy_vdem, 0),
     indeg_democ_power_int_vdem = tidyr::replace_na(indeg_democ_power_int_vdem, 0),
     indeg_women_participation = tidyr::replace_na(indeg_women_participation, 0),
     indeg_autoc_flag = tidyr::replace_na(indeg_autoc_flag, 0), 
     indeg_autoc_power_flag = tidyr::replace_na(indeg_autoc_power_flag, 0), 
     indeg_exec_constraint = tidyr::replace_na(indeg_exec_constraint, 0), 
     indeg_comp_pol_part = tidyr::replace_na(indeg_comp_pol_part, 0), 
     indeg_unli_exec = tidyr::replace_na(indeg_unli_exec, 0), 
     indeg_rest_pol_part = tidyr::replace_na(indeg_rest_pol_part, 0),
     indeg_gov_censorship = tidyr::replace_na(indeg_gov_censorship, 0), 
     indeg_religious_freedom = tidyr::replace_na(indeg_religious_freedom, 0),
     democ_temp = if_else(polity >= 6, 1, 0)
   )
   )  
   
   suppressWarnings(year.data <- mutate(
     year.data,
     indegree = tidyr::replace_na(indegree, 0),
     prop.democ = indeg_democracy,
     prop.democ.mimic = indeg_democ_power_int/indegree,
     prop.women.suffrage.mimic = indeg_women_participation/indegree,
     prop.gov_censorship.mimic = indeg_gov_censorship/indegree,
     prop.religious.freedom.mimic = indeg_religious_freedom/indegree,
     prop.democ.mimic.vdem = indeg_democracy_vdem/indegree,
     prop.autoc = indeg_autoc_flag/indegree,
     prop.autoc.mimic = indeg_autoc_power_flag/indegree,
     prop.exec_constraint = indeg_exec_constraint/indegree,
     prop.comp_pol_part = indeg_comp_pol_part/indegree,
     prop.unli_exec = indeg_unli_exec/indegree,
     prop.rest_pol_part = indeg_rest_pol_part/indegree,
     prop.democ = tidyr::replace_na(prop.democ, 0),
     prop.democ.mimic = tidyr::replace_na(prop.democ.mimic, 0),
     prop.women.suffrage.mimic = tidyr::replace_na(prop.women.suffrage.mimic, 0),
     prop.gov_censorship.mimic = tidyr::replace_na(prop.gov_censorship.mimic, 0),
     prop.religious.freedom.mimic = tidyr::replace_na(prop.religious.freedom.mimic, 0),
     prop.democ.mimic.vdem = tidyr::replace_na(prop.democ.mimic.vdem, 0),
     prop.autoc = tidyr::replace_na(prop.autoc, 0),
     prop.autoc.mimic = tidyr::replace_na(prop.autoc.mimic, 0),
     prop.exec_constraint = tidyr::replace_na(prop.exec_constraint, 0),
     prop.comp_pol_part = tidyr::replace_na(prop.comp_pol_part, 0),
     prop.unli_exec = tidyr::replace_na(prop.unli_exec, 0),
     prop.rest_pol_part = tidyr::replace_na(prop.rest_pol_part, 0),
     prop.democ.mimic.vdem = ifelse(prop.democ.mimic.vdem > 1, 0, prop.democ.mimic.vdem),
     prop.autoc = ifelse(prop.autoc > 1, 0, prop.autoc),
     prop.autoc.mimic = ifelse(prop.autoc.mimic > 1, 0, prop.autoc.mimic),
     prop.exec_constraint = ifelse(prop.exec_constraint > 1, 0, prop.exec_constraint),
     prop.comp_pol_part = ifelse(prop.comp_pol_part > 1, 0, prop.comp_pol_part),
     prop.unli_exec = ifelse(prop.unli_exec > 1, 0, prop.unli_exec),
     prop.rest_pol_part = ifelse(prop.rest_pol_part > 1, 0, prop.rest_pol_part),
     names = to
   ) %>% dplyr::select(
     names, year, prop.democ.mimic, prop.autoc.mimic,
     prop.gov_censorship.mimic, prop.rest_pol_part,
     prop.exec_constraint, v2mecenefm_ord
   ) 
   )
   year.data <- data.frame(year.data)
   
   if (all(is.na(year.data$gle_gdp)))
   {
     year.data[,substr(colnames(year.data), nchar(colnames(year.data)) - 3, nchar(colnames(year.data))) == "_gdp"] <- NA
   }
   
   if (all(is.na(year.data$v2mecenefm_ord)))
   {
     year.data$prop.gov_censorship.mimic <- NA
   }
   
   
   
   polity.diffusion.permutation[[j]] <- year.data
 }
 polity.diffusion.permutation_temp <- do.call(rbind, polity.diffusion.permutation)
 polity.diffusion.permutation_temp <- left_join(
   polity.diffusion.merge,
   polity.diffusion.permutation_temp, 
   by = c(
     "names", "year"
   )
 ) %>% mutate(
   prop.democ.mimic = ifelse(is.na(prop.democ.mimic), 0, prop.democ.mimic), 
   prop.autoc.mimic = ifelse(is.na(prop.autoc.mimic), 0, prop.autoc.mimic), 
   prop.gov_censorship.mimic = ifelse(is.na(prop.gov_censorship.mimic), 0, prop.gov_censorship.mimic), 
   prop.rest_pol_part = ifelse(is.na(prop.rest_pol_part), 0, prop.rest_pol_part), 
   prop.exec_constraint = ifelse(is.na(prop.exec_constraint), 0, prop.exec_constraint)
 )
 permutation_networks[[permute_nets]] <- polity.diffusion.permutation_temp
 cat("\r")
 cat("\rPermute network ties: ", permute_nets)
}
save(permutation_networks, file = "model_output/permutation_networks_emulate.rda")
load("model_output/permutation_networks_emulate.rda")

## Analyze permuted networks
suppressWarnings({
  permutation_values <- list()
  for (i in 1:length(permutation_networks))
  {
    permutation_networks[[i]] <- mutate(
      permutation_networks[[i]],
      prop.democ.mimic = ifelse(is.infinite(prop.democ.mimic), 0, prop.democ.mimic),
      prop.autoc.mimic = ifelse(is.infinite(prop.autoc.mimic), 0, prop.autoc.mimic),
      prop.rest_pol_part = ifelse(is.infinite(prop.rest_pol_part), 0, prop.rest_pol_part),
      prop.gov_censorship.mimic = ifelse(is.infinite(prop.gov_censorship.mimic), 0, prop.gov_censorship.mimic),
      prop.exec_constraint = ifelse(is.infinite(prop.exec_constraint), 0, prop.exec_constraint)
    )
    
    perm.democ         <- glmer(democ ~ prop.democ.mimic + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                              family = binomial(link = "logit"), 
                              control = glmerControl(
                                optimizer = "bobyqa", optCtrl = list(maxfun = 1e10)
                              ), data =  permutation_networks[[i]])
    perm.autoc         <- glmer(autoc ~ prop.autoc.mimic + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                              family = binomial(link = "logit"), 
                              control = glmerControl(
                                optimizer = "bobyqa", optCtrl = list(maxfun = 1e10)
                              ), data =  permutation_networks[[i]])
    prem.rest.pol.part <- glmer(restricted_political_participation ~ prop.rest_pol_part + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                              family = binomial(link = "logit"), 
                              control = glmerControl(
                                optimizer = "bobyqa", optCtrl = list(maxfun = 1e10)
                              ),  data =  permutation_networks[[i]])
    perm.free.media    <- glmer(free_media ~ prop.gov_censorship.mimic + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2) + (1|year) + (1|names), 
                              family = binomial(link = "logit"), 
                              control = glmerControl(
                                optimizer = "bobyqa", optCtrl = list(maxfun = 1e10)
                              ),  data =  permutation_networks[[i]])
    perm.exec.const    <- glmer(executive_constrainted ~ + prop.exec_constraint  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                              family = binomial(link = "logit"), 
                              control = glmerControl(
                                optimizer = "bobyqa", optCtrl = list(maxfun = 1e10)
                              ),  data =  permutation_networks[[i]])
    
    try({perm.democ_value         <- fixef(perm.democ)[2]})
    try({perm.autoc_value         <- fixef(perm.autoc)[2]})
    try({prem.rest.pol.part_value <- fixef(prem.rest.pol.part)[2]})
    try({perm.free.media_value    <- fixef(perm.free.media)[2]})
    try({perm.exec.const_value    <- fixef(perm.exec.const)[2]})
    
    
    try({
      return.item <- data.frame(
        perm_democ = perm.democ_value,
        perm_autoc = perm.autoc_value,
        perm_pol_rest_part = prem.rest.pol.part_value,
        perm_free_media = perm.free.media_value,
        perm_exec_const = perm.exec.const_value
        )})
    
    try({rm(perm.democ_value)})
    try({rm(perm.autoc_value)})
    try({rm(prem.rest.pol.part_value)})
    try({rm(perm.free.media_value)})
    try({rm(perm.exec.const_value)})
    try({row.names(return.item) <- NULL})
    permutation_values[[i]] <- return.item
    cat("\r")
    cat(paste0("\rPercent analyze permuted networks: ", round(i/length(permutation_networks), 2)))
  }
})
permutation_estimates <- do.call(bind_rows, permutation_values)
save(permutation_estimates, file = "model_output/permutation_estimates_emulate.rda")

## Run some robustness checks
rm(list = ls())
range01 <- function(x){(x-min(x, na.rm = T))/(max(x, na.rm = T)-min(x, na.rm = T))}


## Load in data
adjacency_matricies <- readRDS("intl_order/adjacency_matricies_replication.rds")
polity.diffusion <- readRDS("data/polity_diffusion_replication.rds")
min_year <- min(polity.diffusion$year) - 1

polity.diffusion <- mutate(
  polity.diffusion,
  democ = ifelse(polity >= 6, 1, 0),
  autoc = ifelse(polity <= -6, 1, 0),
  v2csgender_ord = as.numeric(as.character(v2csgender_ord)),
  v2csrlgrep_ord = as.numeric(as.character(v2csrlgrep_ord)),
  parcomp = as.numeric(as.character(parcomp)),
  parreg = as.numeric(as.character(parreg)),
  xconst = as.numeric(as.character(xconst)),
  v2mecenefm_ord = as.numeric(as.character(v2mecenefm_ord)),
  women_suffrage = ifelse(v2csgender_ord >= 3, 1, 0),
  religious_freedom = ifelse(v2csrlgrep_ord >= 3, 1, 0),
  free_media = ifelse(v2mecenefm_ord >= 3, 1, 0),
  wdi_gdpcapgr = as.numeric(scale(wdi_gdpcapgr)),
  comptetive_political_participation = ifelse(parcomp == 5, 1, 0),
  unilateral_executive = ifelse(xconst  == 1, 1, 0),
  restricted_political_participation = ifelse(parreg == 4, 1, 0),
  executive_constrainted = ifelse(xconst == 7, 1, 0),
  isolated = ifelse(indeg_power == 0, 1, 0)
)

polity.diffusion$year_num <- as.numeric(scale(polity.diffusion$year))
polity.diffusion$gle_cgdpc <- as.numeric(scale(log(polity.diffusion$gle_cgdpc)))
polity.diffusion$wdi_gdpcapgr <- as.numeric(scale(log(polity.diffusion$wdi_gdpcapgr + -1 * min(polity.diffusion$wdi_gdpcapgr, na.rm = T) + 1)))

set.seed(1)
polity.diffusion$vdem_democ <- ifelse(polity.diffusion$v2x_polyarchy >= 0.5, 1, 0)
re.democ.vdem <- glmer(vdem_democ ~ prop_democ_mimic_vdem + democratic_neighbor_vdem + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                  family = binomial(link = "logit"), 
                  control = glmerControl(optimizer = "bobyqa",
                                         optCtrl = list(maxfun = 5e10)), 
                  data = polity.diffusion)


re.women_suff <- glmer(women_suffrage ~ prop_women_suffrage_mimic + women_participation_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                  family = binomial(link = "logit"), 
                  control = glmerControl(optimizer = "bobyqa",
                                         optCtrl = list(maxfun = 5e10))
                  , 
                  data = polity.diffusion)

re.unil_exec <- glmer(unilateral_executive ~ prop_unli_exec + unli_exec_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                       family = binomial(link = "logit"), 
                       control = glmerControl(optimizer = "bobyqa",
                                              optCtrl = list(maxfun = 5e10))
                       , 
                       data = polity.diffusion)

re.rel.freedom <- glmer(religious_freedom ~ prop_religious_freedom_mimic + religious_freedom_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year, 2) + (1|year) + (1|names), 
                        family = binomial(link = "logit"), 
                        control = glmerControl(optimizer = "bobyqa",
                                               optCtrl = list(maxfun = 5e10))
                        ,  
                        data = polity.diffusion)

point.est.re.democ <- glmer(democ ~ prop_democ_mimic + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                  family = binomial(link = "logit"), 
                  control = glmerControl(
                    optimizer = "bobyqa",
                    optCtrl = list(maxfun = 5e10), 
                    check.scaleX = "ignore",
                    check.conv.hess = "ignore",
                    check.conv.grad = "ignore"), 
                  data = polity.diffusion)

point.est.re.autoc <- glmer(autoc ~ prop_autoc_mimic + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                  family = binomial(link = "logit"), 
                  control = glmerControl(optimizer = "bobyqa",
                                         optCtrl = list(maxfun = 5e10))
                  , 
                  data = polity.diffusion)

point.est.re.rest.part <- glmer(restricted_political_participation ~ prop_rest_pol_part + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                      family = binomial(link = "logit"), 
                      control = glmerControl(optimizer = "bobyqa",
                                             optCtrl = list(maxfun = 5e10))
                      ,  
                      data = polity.diffusion)


point.est.re.free.media <- glmer(free_media ~ prop_gov_censorship_mimic + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2) + (1|year) + (1|names), 
                       family = binomial(link = "logit"), 
                       control = glmerControl(optimizer = "bobyqa",
                                              optCtrl = list(maxfun = 5e10))
                       ,  
                       data = polity.diffusion)


point.est.re.exec.const <- glmer(executive_constrainted ~ prop_exec_constraint + exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                       family = binomial(link = "logit"), 
                       control = glmerControl(
                         optimizer = "bobyqa", optCtrl = list(maxfun = 5e10)
                       ), 
                       data = polity.diffusion)

glm.point.est.re.democ <- glm(democ ~ prop_democ_mimic + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                                family = "binomial", 
                                data = polity.diffusion)


glm.point.est.re.autoc <- glm(autoc ~ prop_autoc_mimic + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                                family = "binomial", 
                                data = polity.diffusion)

glm.point.est.re.rest.part <- glm(restricted_political_participation ~ prop_rest_pol_part + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                                    family = "binomial", 
                                    data = polity.diffusion)

glm.point.est.re.free.media <- glm(free_media ~ prop_gov_censorship_mimic + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2), 
                                     family = "binomial", 
                                     data = polity.diffusion)

glm.point.est.re.exec.const <- glm(executive_constrainted ~ prop_exec_constraint + exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                                 family = "binomial", 
                                 data = polity.diffusion)

biv.democ <- glm(democ ~ prop_democ_mimic, family = "binomial", data = polity.diffusion)
biv.autoc <- glm(autoc ~ prop_autoc_mimic, family = "binomial", data = polity.diffusion)
biv.rest.part <- glm(restricted_political_participation ~ prop_rest_pol_part, family = "binomial", data = polity.diffusion)
biv.free.media <- glm(free_media ~ prop_gov_censorship_mimic, family = "binomial", data = polity.diffusion)
biv.exec.const <- glm(executive_constrainted ~ prop_exec_constraint, family = "binomial", data = polity.diffusion)


lm.point.est.re.democ <- glm(democ ~ prop_democ_mimic + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                              data = polity.diffusion)


lm.point.est.re.autoc <- glm(autoc ~ prop_autoc_mimic + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                              data = polity.diffusion)

lm.point.est.re.rest.part <- glm(restricted_political_participation ~ prop_rest_pol_part + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                                  data = polity.diffusion)

lm.point.est.re.free.media <- glm(free_media ~ prop_gov_censorship_mimic + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2), 
                                   data = polity.diffusion)


lm.point.est.re.exec.const <- glm(executive_constrainted ~ prop_exec_constraint + exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                                 data = polity.diffusion)

save(point.est.re.democ, file = "model_output/point.est.re.democ.rda")
save(point.est.re.autoc, file = "model_output/point.est.re.autoc.rda")
save(point.est.re.rest.part, file = "model_output/point.est.re.rest.part.rda")
save(point.est.re.free.media, file = "model_output/point.est.re.free.media.rda")
save(point.est.re.exec.const, file = "model_output/point.est.re.exec.const.rda")

save(glm.point.est.re.democ, file = "model_output/glm.point.est.re.democ.rda")
save(glm.point.est.re.autoc, file = "model_output/glm.point.est.re.autoc.rda")
save(glm.point.est.re.rest.part, file = "model_output/glm.point.est.re.rest.part.rda")
save(glm.point.est.re.free.media, file = "model_output/glm.point.est.re.free.media.rda")
save(glm.point.est.re.exec.const, file = "model_output/glm.point.est.re.exec.const.rda")

save(lm.point.est.re.democ, file = "model_output/lm.point.est.re.democ.rda")
save(lm.point.est.re.autoc, file = "model_output/lm.point.est.re.autoc.rda")
save(lm.point.est.re.rest.part, file = "model_output/lm.point.est.re.rest.part.rda")
save(lm.point.est.re.free.media, file = "model_output/lm.point.est.re.free.media.rda")
save(lm.point.est.re.exec.const, file = "model_output/lm.point.est.re.exec.const.rda")

save(biv.democ, file = "model_output/biv.democ.rda")
save(biv.autoc, file = "model_output/biv.autoc.rda")
save(biv.rest.part, file = "model_output/biv.rest.part.rda")
save(biv.free.media, file = "model_output/biv.free.media.rda")
save(biv.exec.const, file = "model_output/biv.exec.const.rda")

extracted_diff_data <- data.frame(
  mean = c(
    fixef(re.democ.vdem),
    fixef(re.women_suff),
    fixef(re.unil_exec),
    fixef(re.rel.freedom)
  ),
  names = c(
    names(fixef(re.democ.vdem)),
    names(fixef(re.women_suff)),
    names(fixef(re.unil_exec)),
    names(fixef(re.rel.freedom))
  ),
  se = c(
    summary(re.democ.vdem)$coefficients[,2],
    summary(re.women_suff)$coefficients[,2],
    summary(re.unil_exec)$coefficients[,2],
    summary(re.rel.freedom)$coefficients[,2]
  ),
  dv = c(
    rep("Democracy (V-Dem)", 8),
    rep("Women suffrage", 8),
    rep("Unilateral executive", 8),
    rep("Religious freedom", 8)
  )
) |> filter(
  !grepl("year", names)
)

extracted_diff_data <- mutate(
  extracted_diff_data, 
  names = ifelse(names == "(Intercept)", "Constant", names),
  names = ifelse(names == "civil_war", "Civil War", names),
  names = ifelse(names == "gle_cgdpc", "log of\nGDPC", names),
  names = ifelse(names == "wdi_gdpcapgr", "Annual growth\nin GDPC", names),
  names = ifelse(grepl("neighbor", names), "Pct. neighbors\nw/ attribute", names),
  names = ifelse(grepl("prop", names), "Prop. of in-\ndegree ties", names),
  ub = mean + 1.96 * se,
  lb = mean - 1.96 * se,
  sig = ifelse(lb > 0 | ub < 0, "Statistically\nsignificant", "Not statistically\nsignificant")
)


ggplot(extracted_diff_data, aes(x = mean, xmin = lb, xmax = ub, y = reorder(names, lb), col = sig)) + 
  geom_errorbar(width = 0.2, size = 2) + 
  facet_wrap(~dv, scales = "free_x", ncol = 4) + 
  theme_bw() + 
  labs(y = "", x = "Coefficient estimates")  + 
  geom_vline(xintercept = 0, linetype = "dashed", col = "black", size = 1) + 
  scale_colour_brewer(palette = "Dark2") + 
  theme(axis.text = element_text(size = 15),
        axis.title = element_text(size = 15),
        strip.text = element_text(size = 15),
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size = 15))
ggsave("data_visualization/alt_social_measures.png", width = 16, height = 6, units = "in")  

set.seed(1)
polity.diffusion.merge <- dplyr::select(
  polity.diffusion,
  c(
    democ, autoc, executive_constrainted, free_media , restricted_political_participation,
    democratic_neighbor, autocratic_neighbor, exec_constraint_neighbor, gov_censorship_neighbor, rest_pol_part_neighbor,
    gle_cgdpc, wdi_gdpcapgr, civil_war, year_num, year, names,
    v2mecenefm_ord
  )
)
## Permutation test for CINC
permutate_strength_networks <- list()
X <- 1:500
for (permute_nets in X) 
{
  polity.diffusion.permutation <- list()
  for(j in 1:length(adjacency_matricies)) 
  {
    
    g <- graph_from_adjacency_matrix(adjacency_matricies[[j]], mode = "undirected", weighted = T)
    orig_graph    <- reshape2::melt(adjacency_matricies[[j]])
    orig_graph <- subset(orig_graph, Var1 != Var2)
    orig_graph <- mutate(
      orig_graph, 
      Var1  = as.character(Var1),
      Var2  = as.character(Var2),
      Var1a = if_else(Var1 < Var2, Var1, Var2),
      Var2a = if_else(Var2 < Var1, Var1, Var2),
      Var1 = Var1a, 
      Var2 = Var2a
    ) %>% dplyr::select(
      Var1, Var2, value
    ) %>% distinct(
      Var1, Var2, .keep_all = T
    )
    
    indeg <- data.frame(
      names = names(degree(g)),
      indegree = strength(g, mode = "in")
    )
    
    democ_autoc_data <- filter(
      polity.diffusion, 
      year == j + min_year
    ) %>% 
      distinct(
        names, year, .keep_all = T
      )
    
    power_temp <- distinct(
      democ_autoc_data,
      names, cinc
    )
    power_temp$cinc <- sample(power_temp$cinc)
    weighted_edges <- dplyr::rename(
      orig_graph, 
      to = Var1, 
      from = Var2,
      tie = value
    )
    
    edge.indeg <- left_join(
      weighted_edges,
      power_temp,
      by = c("from" = "names")
    ) %>% left_join(
      power_temp,
      by = c("to" = "names")
    ) %>% mutate(
      democ_flag = if_else(from %in% subset(democ_autoc_data, polity >= 6)$names, 1, 0),
      democ_power_flag = if_else(from %in% subset(democ_autoc_data, polity >= 6)$names & cinc.x > cinc.y, 1, 0),
      democ_flag_vdem = if_else(from %in% subset(democ_autoc_data, v2x_polyarchy >= 0.5)$names, 1, 0),
      democ_power_flag_vdem = if_else(from %in% subset(democ_autoc_data, v2x_polyarchy >= 0.5)$names & cinc.x > cinc.y, 1, 0),
      autoc_flag = if_else(from %in% subset(democ_autoc_data, polity <= -6)$names, 1, 0),
      autoc_power_flag = if_else(from %in% subset(democ_autoc_data, polity <= -6)$names & cinc.x > cinc.y, 1, 0),
      exec_constraint = if_else(from %in% subset(democ_autoc_data,  xconst  == 7)$names & cinc.x > cinc.y, 1, 0),
      comp_pol_part = if_else(from %in% subset(democ_autoc_data,  parcomp  == 5)$names & cinc.x > cinc.y, 1, 0),
      unli_exec = if_else(from %in% subset(democ_autoc_data,  xconst  == 1)$names & cinc.x > cinc.y, 1, 0),
      rest_pol_part = if_else(from %in% subset(democ_autoc_data,  parreg  == 4)$names  & cinc.x > cinc.y, 1, 0),
      women_participation_flag = if_else(from %in% subset(democ_autoc_data, v2csgender_ord >= 3)$names & cinc.x > cinc.y, 1, 0),
      gov_censorship_flag = if_else(from %in% subset(democ_autoc_data, v2mecenefm_ord >= 3)$names & cinc.x > cinc.y, 1, 0),
      religious_freedom_flag = if_else(from %in% subset(democ_autoc_data, v2csrlgrep_ord >= 3)$names & cinc.x > cinc.y, 1, 0)
    ) %>% plyr::ddply(
      ~to,
      summarize, 
      indeg_democracy = sum(democ_flag * tie),
      indeg_democ_power_int = sum(democ_power_flag  * tie),
      indeg_democracy_vdem = sum(democ_flag_vdem  * tie),
      indeg_democ_power_int_vdem = sum(democ_power_flag_vdem * tie),
      indeg_women_participation = sum(women_participation_flag * tie),
      indeg_autoc_flag = sum(autoc_flag * tie),
      indeg_autoc_power_flag = sum(autoc_power_flag * tie),
      indeg_exec_constraint = sum(exec_constraint * tie),
      indeg_comp_pol_part = sum(comp_pol_part * tie),
      indeg_unli_exec = sum(unli_exec * tie),
      indeg_rest_pol_part = sum(rest_pol_part * tie),
      indeg_gov_censorship = sum(gov_censorship_flag * tie),
      indeg_religious_freedom = sum(religious_freedom_flag * tie)
    ) %>% left_join(
      indeg, 
      by = c(
        "to" = "names"
      )
    )
    edge.indeg <- edge.indeg[!duplicated(edge.indeg),]
    democ_autoc_data <- data.frame(democ_autoc_data)
    democ_autoc_data <- democ_autoc_data[,!(colnames(democ_autoc_data) %in% c(colnames(edge.indeg), colnames(indeg))[c(colnames(edge.indeg), colnames(indeg)) != "names"])]
    suppressWarnings(year.data <- left_join(
      edge.indeg, 
      democ_autoc_data,
      by = c("to" = "names")
    )  %>% mutate(
      indeg_democracy = tidyr::replace_na(indeg_democracy, 0),
      indeg_democ_power_int = tidyr::replace_na(indeg_democ_power_int, 0),
      indeg_democracy_vdem = tidyr::replace_na(indeg_democracy_vdem, 0),
      indeg_democ_power_int_vdem = tidyr::replace_na(indeg_democ_power_int_vdem, 0),
      indeg_women_participation = tidyr::replace_na(indeg_women_participation, 0),
      indeg_autoc_flag = tidyr::replace_na(indeg_autoc_flag, 0), 
      indeg_autoc_power_flag = tidyr::replace_na(indeg_autoc_power_flag, 0), 
      indeg_exec_constraint = tidyr::replace_na(indeg_exec_constraint, 0), 
      indeg_comp_pol_part = tidyr::replace_na(indeg_comp_pol_part, 0), 
      indeg_unli_exec = tidyr::replace_na(indeg_unli_exec, 0), 
      indeg_rest_pol_part = tidyr::replace_na(indeg_rest_pol_part, 0),
      indeg_gov_censorship = tidyr::replace_na(indeg_gov_censorship, 0), 
      indeg_religious_freedom = tidyr::replace_na(indeg_religious_freedom, 0),
      democ_temp = if_else(polity >= 6, 1, 0)
    )
    )  
    
    suppressWarnings(year.data <- mutate(
      year.data,
      indegree = tidyr::replace_na(indegree, 0),
      prop.democ = indeg_democracy,
      prop.democ.mimic = indeg_democ_power_int/indegree,
      prop.women.suffrage.mimic = indeg_women_participation/indegree,
      prop.gov_censorship.mimic = indeg_gov_censorship/indegree,
      prop.religious.freedom.mimic = indeg_religious_freedom/indegree,
      prop.democ.mimic.vdem = indeg_democracy_vdem/indegree,
      prop.autoc = indeg_autoc_flag/indegree,
      prop.autoc.mimic = indeg_autoc_power_flag/indegree,
      prop.exec_constraint = indeg_exec_constraint/indegree,
      prop.comp_pol_part = indeg_comp_pol_part/indegree,
      prop.unli_exec = indeg_unli_exec/indegree,
      prop.rest_pol_part = indeg_rest_pol_part/indegree,
      prop.democ = tidyr::replace_na(prop.democ, 0),
      prop.democ.mimic = tidyr::replace_na(prop.democ.mimic, 0),
      prop.women.suffrage.mimic = tidyr::replace_na(prop.women.suffrage.mimic, 0),
      prop.gov_censorship.mimic = tidyr::replace_na(prop.gov_censorship.mimic, 0),
      prop.religious.freedom.mimic = tidyr::replace_na(prop.religious.freedom.mimic, 0),
      prop.democ.mimic.vdem = tidyr::replace_na(prop.democ.mimic.vdem, 0),
      prop.democ = tidyr::replace_na(prop.democ, 0),
      prop.democ.mimic = tidyr::replace_na(prop.democ.mimic, 0),
      prop.autoc = tidyr::replace_na(prop.autoc, 0),
      prop.autoc.mimic = tidyr::replace_na(prop.autoc.mimic, 0),
      prop.exec_constraint = tidyr::replace_na(prop.exec_constraint, 0),
      prop.comp_pol_part = tidyr::replace_na(prop.comp_pol_part, 0),
      prop.unli_exec = tidyr::replace_na(prop.unli_exec, 0),
      prop.rest_pol_part = tidyr::replace_na(prop.rest_pol_part, 0),
      prop.democ.mimic.vdem = ifelse(prop.democ.mimic.vdem > 1, 0, prop.democ.mimic.vdem),
      prop.autoc = ifelse(prop.autoc > 1, 0, prop.autoc),
      prop.autoc.mimic = ifelse(prop.autoc.mimic > 1, 0, prop.autoc.mimic),
      prop.exec_constraint = ifelse(prop.exec_constraint > 1, 0, prop.exec_constraint),
      prop.comp_pol_part = ifelse(prop.comp_pol_part > 1, 0, prop.comp_pol_part),
      prop.unli_exec = ifelse(prop.unli_exec > 1, 0, prop.unli_exec),
      prop.rest_pol_part = ifelse(prop.rest_pol_part > 1, 0, prop.rest_pol_part),
      names = to
    )  %>% dplyr::select(
      names, year, prop.democ.mimic, prop.autoc.mimic,
      prop.gov_censorship.mimic, prop.rest_pol_part,
      prop.exec_constraint, v2mecenefm_ord
    ) 
    )
    year.data <- data.frame(year.data)
    
    if (all(is.na(year.data$gle_gdp)))
    {
      year.data[,substr(colnames(year.data), nchar(colnames(year.data)) - 3, nchar(colnames(year.data))) == "_gdp"] <- NA
    }
    
    if (all(is.na(year.data$v2mecenefm_ord)))
    {
      year.data$prop.gov_censorship.mimic <- NA
    }
    
    if (all(is.na(year.data$v2csrlgrep_ord)))
    {
      year.data$prop.religious.freedom.mimic <- NA
    }
    
    polity.diffusion.permutation[[j]] <- year.data
  }
  polity.diffusion.permutation_temp <- do.call(rbind, polity.diffusion.permutation)
  polity.diffusion.permutation_temp <- left_join(
    polity.diffusion,
    polity.diffusion.permutation_temp, 
    by = c(
      "names", "year"
    )
  ) %>% mutate(
    prop.democ.mimic = ifelse(is.na(prop.democ.mimic), 0, prop.democ.mimic), 
    prop.autoc.mimic = ifelse(is.na(prop.autoc.mimic), 0, prop.autoc.mimic), 
    prop.gov_censorship.mimic = ifelse(is.na(prop.gov_censorship.mimic), 0, prop.gov_censorship.mimic), 
    prop.rest_pol_part = ifelse(is.na(prop.rest_pol_part), 0, prop.rest_pol_part), 
    prop.exec_constraint = ifelse(is.na(prop.exec_constraint), 0, prop.exec_constraint)
  )
  
  permutate_strength_networks[[permute_nets]] <- polity.diffusion.permutation_temp
  cat("\rPermute CINC networks: ", permute_nets)
}
save(permutate_strength_networks, file = "model_output/permutate_strength_networks.rda")

load("model_output/permutate_strength_networks.rda")

suppressWarnings({
  permutate_strength_values <- list()
  for (i in 1:length(permutate_strength_networks))
  {
    
    permutate_strength_networks[[i]] <- mutate(
      permutate_strength_networks[[i]],
        prop.democ.mimic = ifelse(is.infinite(prop.democ.mimic), 0, prop.democ.mimic),
        prop.autoc.mimic = ifelse(is.infinite(prop.autoc.mimic), 0, prop.autoc.mimic),
        prop.rest_pol_part = ifelse(is.infinite(prop.rest_pol_part), 0, prop.rest_pol_part),
        prop.gov_censorship.mimic = ifelse(is.infinite(prop.gov_censorship.mimic), 0, prop.gov_censorship.mimic),
        prop.exec_constraint = ifelse(is.infinite(prop.exec_constraint), 0, prop.exec_constraint)
    )
    perm.democ         <- glmer(democ ~ prop.democ.mimic + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                                family = binomial(link = "logit"), 
                                control = glmerControl(
                                  optimizer = "bobyqa", optCtrl = list(maxfun = 1e10)
                                ), data =  permutate_strength_networks[[i]])
    perm.autoc         <- glmer(autoc ~ prop.autoc.mimic + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                                family = binomial(link = "logit"), 
                                control = glmerControl(
                                  optimizer = "bobyqa", optCtrl = list(maxfun = 1e10)
                                ), data =  permutate_strength_networks[[i]])
    prem.rest.pol.part <- glmer(restricted_political_participation ~ prop.rest_pol_part + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                                family = binomial(link = "logit"), 
                                control = glmerControl(
                                  optimizer = "bobyqa", optCtrl = list(maxfun = 1e10)
                                ),  data =  permutate_strength_networks[[i]])
    perm.free.media    <- glmer(free_media ~ prop.gov_censorship.mimic + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2) + (1|year) + (1|names), 
                                family = binomial(link = "logit"), 
                                control = glmerControl(
                                  optimizer = "bobyqa", optCtrl = list(maxfun = 1e10)
                                ),  data =  permutate_strength_networks[[i]])
    perm.exec.const    <- glmer(executive_constrainted ~ + prop.exec_constraint  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                                family = binomial(link = "logit"), 
                                control = glmerControl(
                                  optimizer = "bobyqa", optCtrl = list(maxfun = 1e10)
                                ),  data =  permutate_strength_networks[[i]])
    
    try({perm.democ_value         <- fixef(perm.democ)[2]})
    try({perm.autoc_value         <- fixef(perm.autoc)[2]})
    try({prem.rest.pol.part_value <- fixef(prem.rest.pol.part)[2]})
    try({perm.free.media_value    <- fixef(perm.free.media)[2]})
    try({perm.exec.const_value    <- fixef(perm.exec.const)[2]})
    
    try({
      return.item <- data.frame(
        perm_democ = perm.democ_value,
        perm_autoc = perm.autoc_value,
        perm_pol_rest_part = prem.rest.pol.part_value,
        perm_free_media = perm.free.media_value,
        perm_exec_const = perm.exec.const_value
      )})
    
    try({rm(perm.democ_value)})
    try({rm(perm.autoc_value)})
    try({rm(prem.rest.pol.part_value)})
    try({rm(perm.free.media_value)})
    try({rm(perm.exec.const_value)})
    
    try({row.names(return.item) <- NULL})
    permutate_strength_values[[i]] <- return.item
    cat("\r")
    cat(paste0("\rPercent analyze permuted networks: ", round(i/length(permutate_strength_networks), 2)))
  }
})

save(permutate_strength_values, file = "model_output/permutate_strength_values.rda")
permutation_strength_estimates <- do.call(bind_rows, permutate_strength_values)
save(permutation_strength_estimates, file = "model_output/permutate_strength_estimates.rda")

rm(list = ls())
range01 <- function(x){(x-min(x, na.rm = T))/(max(x, na.rm = T)-min(x, na.rm = T))}

## Load in data
adjacency_matricies <- readRDS("intl_order/adjacency_matricies_replication.rds")
polity.diffusion <- readRDS("data/polity_diffusion_replication.rds")
min_year <- min(polity.diffusion$year) - 1
polity.diffusion$year_num <- as.numeric(scale(polity.diffusion$year))
polity.diffusion$gle_cgdpc <- as.numeric(scale(log(polity.diffusion$gle_cgdpc)))
polity.diffusion$wdi_gdpcapgr <- as.numeric(scale(log(polity.diffusion$wdi_gdpcapgr + -1 * min(polity.diffusion$wdi_gdpcapgr, na.rm = T) + 1)))

## Create variables
polity.diffusion <- mutate(
  polity.diffusion,
  democ = ifelse(polity >= 6, 1, 0),
  autoc = ifelse(polity <= -6, 1, 0),
  v2csgender_ord = as.numeric(as.character(v2csgender_ord)),
  v2csrlgrep_ord = as.numeric(as.character(v2csrlgrep_ord)),
  parcomp = as.numeric(as.character(parcomp)),
  parreg = as.numeric(as.character(parreg)),
  xconst = as.numeric(as.character(xconst)),
  v2mecenefm_ord = as.numeric(as.character(v2mecenefm_ord)),
  women_suffrage = ifelse(v2csgender_ord >= 3, 1, 0),
  religious_freedom = ifelse(v2csrlgrep_ord >= 3, 1, 0),
  free_media = ifelse(v2mecenefm_ord >= 3, 1, 0),
  wdi_gdpcapgr = as.numeric(scale(wdi_gdpcapgr)),
  comptetive_political_participation = ifelse(parcomp == 5, 1, 0),
  unilateral_executive = ifelse(xconst  == 1, 1, 0),
  restricted_political_participation = ifelse(parreg == 4, 1, 0),
  executive_constrainted = ifelse(xconst == 7, 1, 0)
)

## Institutional diffusion of democracy
suppressWarnings({
  gdp.re.democ <- glmer(democ ~ prop_democ_mimic_gdp + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                    family = binomial(link = "logit"), 
                    control = glmerControl(
                      optimizer = "bobyqa",
                      optCtrl = list(maxfun = 5e10), 
                      check.scaleX = "ignore",
                      check.conv.hess = "ignore",
                      check.conv.grad = "ignore"), 
                    data = polity.diffusion)
  
  
  
  gdp.glm.democ <- glm(democ ~ prop_democ_mimic_gdp + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                   family = "binomial", 
                   data = polity.diffusion)
  
  gdp.lm.democ <- glm(democ ~ prop_democ_mimic_gdp + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                       data = polity.diffusion)
  
  
  ## Institutional diffusion of autocracy
  gdp.re.autoc <- glmer(autoc ~ prop_autoc_mimic_gdp + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                    family = binomial(link = "logit"), 
                    control = glmerControl(optimizer = "bobyqa",
                                           optCtrl = list(maxfun = 5e10))
                    , 
                    data = polity.diffusion)
  
  
  gdp.glm.autoc <- glm(autoc ~ prop_autoc_mimic_gdp + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                   family = "binomial", 
                   data = polity.diffusion)
  
  gdp.lm.autoc <- glm(autoc ~ prop_autoc_mimic_gdp + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                       data = polity.diffusion)
  
  ## Institutional diffusion of political participation restrictions
  gdp.re.rest.part <- glmer(restricted_political_participation ~ prop_rest_pol_part_gdp + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                        family = binomial(link = "logit"), 
                        control = glmerControl(optimizer = "bobyqa",
                                               optCtrl = list(maxfun = 5e10))
                        ,  
                        data = polity.diffusion)
  
  
  gdp.glm.rest.part <- glm(restricted_political_participation ~ prop_rest_pol_part_gdp + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                       family = "binomial",  
                       data = polity.diffusion)
  
  gdp.lm.rest.part <- glm(restricted_political_participation ~ prop_rest_pol_part_gdp + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                           data = polity.diffusion)
  
  ## Institutional diffusion of free media
  gdp.re.free.media <- glmer(free_media ~ prop_gov_censorship_mimic_gpd + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2) + (1|year) + (1|names), 
                         family = binomial(link = "logit"), 
                         control = glmerControl(optimizer = "bobyqa",
                                                optCtrl = list(maxfun = 5e10))
                         ,  
                         data = polity.diffusion)
  
  
  gdp.glm.free.media <- glm(free_media ~ prop_gov_censorship_mimic_gpd + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2),
                        family = "binomial",
                        data = polity.diffusion)
  
  gdp.lm.free.media <- glm(free_media ~ prop_gov_censorship_mimic_gpd + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2),
                            data = polity.diffusion)
  
  ## Institutional diffusion of executive constraints
  gdp.re.exec.const <- glmer(executive_constrainted ~ prop_exec_constraint_gdp + exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                         family = binomial(link = "logit"), 
                         control = glmerControl(
                           optimizer = "bobyqa", optCtrl = list(maxfun = 5e10)
                         ), 
                         data = polity.diffusion)
  
  gdp.glm.exec.const <- glm(executive_constrainted ~ prop_exec_constraint_gdp + exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                        family = "binomial", 
                        data = polity.diffusion)
  
  gdp.lm.exec.const <- glm(executive_constrainted ~ prop_exec_constraint_gdp + exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                            data = polity.diffusion)
  
})

save(gdp.re.democ, file = "model_output/gdp.re.democ.rda")
save(gdp.re.autoc, file = "model_output/gdp.re.autoc.rda")
save(gdp.re.rest.part, file = "model_output/gdp.re.rest.part.rda")
save(gdp.re.free.media, file = "model_output/gdp.re.free.media.rda")
save(gdp.re.exec.const, file = "model_output/gdp.re.exec.const.rda")

save(gdp.glm.democ, file = "model_output/gdp.glm.democ.rda")
save(gdp.glm.autoc, file = "model_output/gdp.glm.autoc.rda")
save(gdp.glm.rest.part, file = "model_output/gdp.glm.rest.part.rda")
save(gdp.glm.free.media, file = "model_output/gdp.glm.free.media.rda")
save(gdp.glm.exec.const, file = "model_output/gdp.glm.exec.const.rda")

save(gdp.lm.democ, file = "model_output/gdp.lm.democ.rda")
save(gdp.lm.autoc, file = "model_output/gdp.lm.autoc.rda")
save(gdp.lm.rest.part, file = "model_output/gdp.lm.rest.part.rda")
save(gdp.lm.free.media, file = "model_output/gdp.lm.free.media.rda")
save(gdp.lm.exec.const, file = "model_output/gdp.lm.exec.const.rda")

polity.diffusion$democ                              <- factor(polity.diffusion$democ)
polity.diffusion$autoc                              <- factor(polity.diffusion$autoc)
polity.diffusion$executive_constrainted             <- factor(polity.diffusion$executive_constrainted)
polity.diffusion$free_media                         <- factor(polity.diffusion$free_media)
polity.diffusion$restricted_political_participation <- factor(polity.diffusion$restricted_political_participation)


set.seed(1)
lambda_grid <- seq(0, 1, 0.25)
alpha_grid <- seq(1)
srchGrid <- expand.grid(alpha = alpha_grid, lambda = lambda_grid)

polity.diffusion.democ.enet <- dplyr::select(
  polity.diffusion,
  c(
    democ, prop_democ_mimic, democratic_neighbor, 
    gle_cgdpc, wdi_gdpcapgr, civil_war, year_num
  )
) 
polity.diffusion.democ.enet <- polity.diffusion.democ.enet[complete.cases(polity.diffusion.democ.enet),]
democ_elnet = train(
  democ ~ prop_democ_mimic + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
  data = polity.diffusion.democ.enet,
  method = "glmnet",
  tuneGrid = srchGrid,
  allowParallel = TRUE
)

polity.diffusion.exec.const.enet <- dplyr::select(
  polity.diffusion,
  c(
    executive_constrainted, prop_exec_constraint, exec_constraint_neighbor, 
    gle_cgdpc, wdi_gdpcapgr, civil_war, year_num
  )
) 
polity.diffusion.exec.const.enet <- polity.diffusion.exec.const.enet[complete.cases(polity.diffusion.exec.const.enet),]
exec_const_elnet = train(
  executive_constrainted ~ prop_exec_constraint + exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
  data = polity.diffusion.exec.const.enet,
  method = "glmnet",
  tuneGrid = srchGrid,
  allowParallel = TRUE
)

polity.diffusion.autoc.enet <- dplyr::select(
  polity.diffusion,
  c(
    autoc, prop_autoc_mimic, autocratic_neighbor, 
    gle_cgdpc, wdi_gdpcapgr, civil_war, year_num
  )
) 
polity.diffusion.autoc.enet <- polity.diffusion.autoc.enet[complete.cases(polity.diffusion.autoc.enet),]
autoc_elnet = train(
  autoc ~ prop_autoc_mimic + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
  data = polity.diffusion.autoc.enet,
  method = "glmnet",
  tuneGrid = srchGrid,
  allowParallel = TRUE
)

polity.diffusion.free.media.enet <- dplyr::select(
  polity.diffusion,
  c(
    free_media, prop_gov_censorship_mimic, gov_censorship_neighbor,
    gle_cgdpc, wdi_gdpcapgr, civil_war, year_num
  )
) 
polity.diffusion.free.media.enet <- polity.diffusion.free.media.enet[complete.cases(polity.diffusion.free.media.enet),]
free_media_elnet = train(
  free_media ~ prop_gov_censorship_mimic + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2), 
  data = polity.diffusion.free.media.enet,
  method = "glmnet",
  tuneGrid = srchGrid,
  allowParallel = TRUE
)

polity.diffusion.rest.part.enet <- dplyr::select(
  polity.diffusion,
  c(
    restricted_political_participation, prop_rest_pol_part, rest_pol_part_neighbor,
    gle_cgdpc, wdi_gdpcapgr, civil_war, year_num
  )
) 
polity.diffusion.rest.part.enet <- polity.diffusion.rest.part.enet[complete.cases(polity.diffusion.rest.part.enet),]
rest_part_elnet = train(
  restricted_political_participation ~ prop_rest_pol_part + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
  data = polity.diffusion.rest.part.enet,
  method = "glmnet",
  tuneGrid = srchGrid,
  allowParallel = TRUE
)

save(exec_const_elnet, file = "model_output/exec_const_elnet.rda")
save(free_media_elnet, file = "model_output/free_media_elnet.rda")
save(autoc_elnet, file = "model_output/autoc_elnet.rda")
save(democ_elnet, file = "model_output/democ_elnet.rda")
save(rest_part_elnet, file = "model_output/rest_part_elnet.rda")


## Last replication unweighted
polity.unimputed <- readRDS("data/unimputed_values_replication.rds")
unweighted_data   <- readRDS("intl_order/unweighted_data_replication.rds")

# Contiguity
contiguity_temp <- distinct(
  unweighted_data,
  cowc1, cowc2, year, contig
)

contiguity <- distinct(
  unweighted_data,
  cowc1, cowc2, year, contig
) %>% mutate(
  ccode1a = cowc2, 
  ccode2a = cowc1, 
  cowc1 = ccode1a,
  cowc2 = ccode2a
) %>% dplyr::select(
  cowc1, cowc2, year, contig
) %>% bind_rows(
  contiguity_temp
)

polity.diffusion.unimputed <- list()
for (j in 1:length(adjacency_matricies))
{
  
  ## Extract network
  g <- graph_from_adjacency_matrix(adjacency_matricies[[j]], mode = "directed", weighted = T)
  
  ## Weighted tie strength
  indeg <- data.frame(
    names = names(degree(g)),
    indegree = strength(g, mode = "in")
  )
  
  ## Turn into edge list
  weighted_edges <- reshape2::melt(adjacency_matricies[[j]])
  
  weighted_edges <- filter(
    weighted_edges, 
    Var1 != Var2 ## Drop self loops
  ) %>% dplyr::rename(
    from = Var1, ## Sender
    to   = Var2, ## Reciever
    tie = value  ## Edge strength
  )
  
  ## Subset data
  democ_autoc_data <- filter(
    polity.unimputed, 
    year == j + min_year
  ) %>% 
    distinct(
      names, year, .keep_all = T
    )
  
  ## Make sure it is unique
  power_temp <- distinct(
    democ_autoc_data,
    names, cinc
  )
  
  ## Merge in power
  edge.indeg <- left_join(
    weighted_edges,
    power_temp,
    by = c("from" = "names")
  ) %>% left_join(
    power_temp,
    by = c("to" = "names")
  )  %>% mutate(
    power_flag = ifelse(cinc.x > cinc.y, 1, 0),
    democ_flag = ifelse(from %in% subset(democ_autoc_data, polity >= 6)$names, 1, 0),
    democ_power_flag = ifelse(from %in% subset(democ_autoc_data, polity >= 6)$names & cinc.x > cinc.y, 1, 0),
    democ_power_flag_vdem = ifelse(from %in% subset(democ_autoc_data, v2x_polyarchy >= 0.5)$names & cinc.x > cinc.y, 1, 0),
    autoc_flag = ifelse(from %in% subset(democ_autoc_data, polity <= -6)$names, 1, 0),
    autoc_power_flag = ifelse(from %in% subset(democ_autoc_data, polity <= -6)$names & cinc.x > cinc.y, 1, 0),
    exec_constraint = ifelse(from %in% subset(democ_autoc_data,  xconst  == 7)$names & cinc.x > cinc.y, 1, 0),
    rest_pol_part = ifelse(from %in% subset(democ_autoc_data,  parreg  == 4)$names & cinc.x > cinc.y, 1, 0),
    gov_censorship_flag = ifelse(from %in% subset(democ_autoc_data, v2mecenefm_ord >= 3)$names & cinc.x > cinc.y, 1, 0)
  ) %>% plyr::ddply(
    ~to,
    summarize, 
    indeg_power = sum(power_flag * tie),
    indeg_democracy = sum(democ_flag * tie),
    indeg_democ_power_int = sum(democ_power_flag  * tie),
    indeg_autoc_flag = sum(autoc_flag * tie),
    indeg_autoc_power_flag = sum(autoc_power_flag * tie),
    indeg_exec_constraint = sum(exec_constraint * tie),
    indeg_rest_pol_part = sum(rest_pol_part * tie),
    indeg_gov_censorship = sum(gov_censorship_flag * tie)
  ) 
  
  edge.indeg <- edge.indeg[!duplicated(edge.indeg),]
  
  ## Find neighbors with attribute
  neighbor.status <- filter(
    contiguity,
    year == j + min_year,
    contig == 1
  ) %>% 
    rename(
      from = cowc1,
      to = cowc2
    )
  
  ## Create neighbor with attribute variable
  neighbor.status <- mutate(
    neighbor.status,
    n = 1,
    democratic_neighbor = ifelse(from %in% subset(democ_autoc_data,  polity >= 6)$names, 1, 0),
    democratic_neighbor_vdem = ifelse(from %in% subset(democ_autoc_data,  v2x_polyarchy >= 0.5)$names, 1, 0),
    autocratic_neighbor = ifelse(from %in% subset(democ_autoc_data,  polity <= -6)$names, 1, 0),
    exec_constraint_neighbor = ifelse(from %in% subset(democ_autoc_data,  xconst  == 7)$names, 1, 0),
    comp_pol_part_neighbor = ifelse(from %in% subset(democ_autoc_data,  parcomp  == 5)$names, 1, 0),
    unli_exec_neighbor = ifelse(from %in% subset(democ_autoc_data,  xconst  == 1)$names, 1, 0),
    rest_pol_part_neighbor = ifelse(from %in% subset(democ_autoc_data,  parreg  == 4)$names, 1, 0),
    women_participation_neighbor = ifelse(from %in% subset(democ_autoc_data, v2csgender_ord >= 3)$names, 1, 0),
    gov_censorship_neighbor = ifelse(from %in% subset(democ_autoc_data, v2mecenefm_ord >= 3)$names, 1, 0),
    religious_freedom_neighbor = ifelse(from %in% subset(democ_autoc_data, v2csrlgrep_ord >= 3)$names, 1, 0)
  ) %>% 
    plyr::ddply(
      ~to, 
      summarize, 
      democratic_neighbor = sum(democratic_neighbor)/sum(n),
      autocratic_neighbor = sum(autocratic_neighbor)/sum(n),
      democratic_neighbor_vdem = sum(democratic_neighbor_vdem)/sum(n),
      exec_constraint_neighbor = sum(exec_constraint_neighbor)/sum(n),
      comp_pol_part_neighbor = sum(comp_pol_part_neighbor)/sum(n),
      unli_exec_neighbor = sum(unli_exec_neighbor)/sum(n),
      rest_pol_part_neighbor = sum(rest_pol_part_neighbor)/sum(n),
      women_participation_neighbor = sum(women_participation_neighbor)/sum(n),
      gov_censorship_neighbor = sum(gov_censorship_neighbor)/sum(n),
      religious_freedom_neighbor = sum(religious_freedom_neighbor)/sum(n),
      democratic_neighbor = tidyr::replace_na(democratic_neighbor, 0),
      autocratic_neighbor = tidyr::replace_na(autocratic_neighbor, 0),
      democratic_neighbor_vdem = tidyr::replace_na(democratic_neighbor_vdem, 0),
      exec_constraint_neighbor = tidyr::replace_na(exec_constraint_neighbor, 0),
      comp_pol_part_neighbor = tidyr::replace_na(comp_pol_part_neighbor, 0),
      unli_exec_neighbor = tidyr::replace_na(unli_exec_neighbor, 0),
      rest_pol_part_neighbor = tidyr::replace_na(rest_pol_part_neighbor, 0),
      women_participation_neighbor = tidyr::replace_na(women_participation_neighbor, 0)
    )
  
  ## Join the data together
  suppressWarnings(year.data <- left_join(
    indeg,
    edge.indeg, 
    by = c(
      "names" = "to"
    )                                      
  ) %>% left_join(
    neighbor.status,
    by = c(
      "names" = "to"
    )
  ) %>% left_join(
    democ_autoc_data,
    by = c("names")
  )  %>% mutate(
    indeg_democracy = tidyr::replace_na(indeg_democracy, 0),
    indeg_democ_power_int = tidyr::replace_na(indeg_democ_power_int, 0),
    indeg_autoc_flag = tidyr::replace_na(indeg_autoc_flag, 0), 
    indeg_autoc_power_flag = tidyr::replace_na(indeg_autoc_power_flag, 0), 
    indeg_exec_constraint = tidyr::replace_na(indeg_exec_constraint, 0), 
    indeg_rest_pol_part = tidyr::replace_na(indeg_rest_pol_part, 0),
    indeg_gov_censorship = tidyr::replace_na(indeg_gov_censorship, 0),
    democ_temp = ifelse(polity >= 6, 1, 0),
    democratic_neighbor = tidyr::replace_na(democratic_neighbor, 0),
    autocratic_neighbor = tidyr::replace_na(autocratic_neighbor, 0),
    democratic_neighbor_vdem = tidyr::replace_na(democratic_neighbor_vdem, 0),
    exec_constraint_neighbor = tidyr::replace_na(exec_constraint_neighbor, 0),
    comp_pol_part_neighbor = tidyr::replace_na(comp_pol_part_neighbor, 0),
    unli_exec_neighbor = tidyr::replace_na(unli_exec_neighbor, 0),
    rest_pol_part_neighbor = tidyr::replace_na(rest_pol_part_neighbor, 0)
  )  )
  
  ## Fix missing data
  suppressWarnings(year.data <- mutate(
    year.data,
    indegree = tidyr::replace_na(indegree, 0),
    prop.democ = indeg_democracy/indegree,
    prop.democ.mimic = indeg_democ_power_int/indegree,
    prop.gov_censorship.mimic = indeg_gov_censorship/indegree,
    prop.autoc = indeg_autoc_flag/indegree,
    prop.autoc.mimic = indeg_autoc_power_flag/indegree,
    prop.exec_constraint = indeg_exec_constraint/indegree,
    prop.rest_pol_part = indeg_rest_pol_part/indegree,
    prop.democ = tidyr::replace_na(prop.democ, 0),
    prop.democ.mimic = tidyr::replace_na(prop.democ.mimic, 0),
    prop.gov_censorship.mimic = tidyr::replace_na(prop.gov_censorship.mimic, 0),
    prop.autoc = tidyr::replace_na(prop.autoc, 0),
    prop.autoc.mimic = tidyr::replace_na(prop.autoc.mimic, 0),
    prop.exec_constraint = tidyr::replace_na(prop.exec_constraint, 0),
    prop.rest_pol_part = tidyr::replace_na(prop.rest_pol_part, 0)
  ) 
  )
  colnames(year.data) <- gsub("\\.", "_", colnames(year.data))
  colnames(year.data) <- trimws(colnames(year.data))
  
  if (all(is.na(year.data$gle_gdp)))
  {
    year.data[,substr(colnames(year.data), nchar(colnames(year.data)) - 3, nchar(colnames(year.data))) == "_gdp"] <- NA
  }
  
  if (all(is.na(year.data$v2mecenefm_ord)))
  {
    year.data$prop_gov_censorship_mimic <- NA
  }
  
  if (all(is.na(year.data$v2csrlgrep_ord)))
  {
    year.data$prop_religious_freedom_mimic <- NA
  }
  
  if (all(is.na(year.data$v2x_polyarchy)))
  {
    year.data$prop_democ_mimic_vdem <- NA
  }
  
  if (all(is.na(year.data$v2csgender_ord)))
  {
    year.data$prop_women_suffrage_mimic <- NA
  }
  
  
  polity.diffusion.unimputed[[j]] <- year.data
}

## Bring together data
polity.diffusion.unimputed <- do.call(bind_rows, polity.diffusion.unimputed)


## Create variables
polity.diffusion.unimputed <- mutate(
  polity.diffusion.unimputed,
  democ = ifelse(polity >= 6, 1, 0),
  autoc = ifelse(polity <= -6, 1, 0),
  v2csgender_ord = as.numeric(as.character(v2csgender_ord)),
  v2csrlgrep_ord = as.numeric(as.character(v2csrlgrep_ord)),
  parcomp = as.numeric(as.character(parcomp)),
  parreg = as.numeric(as.character(parreg)),
  xconst = as.numeric(as.character(xconst)),
  v2mecenefm_ord = as.numeric(as.character(v2mecenefm_ord)),
  free_media = ifelse(v2mecenefm_ord >= 3, 1, 0),
  restricted_political_participation = ifelse(parreg == 4, 1, 0),
  executive_constrainted = ifelse(xconst == 7, 1, 0)
)

polity.diffusion.unimputed <- dplyr::select(
  polity.diffusion,
  c(
    names, year, gle_cgdpc, wdi_gdpcapgr, civil_war, year_num
  )
) %>% distinct(
) %>% right_join(
  polity.diffusion.unimputed
)

suppressWarnings({
  unimputed.re.democ <- glmer(democ ~ prop_democ_mimic + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                              family = binomial(link = "logit"), 
                              control = glmerControl(
                                optimizer = "bobyqa",
                                optCtrl = list(maxfun = 5e10), 
                                check.scaleX = "ignore",
                                check.conv.hess = "ignore",
                                check.conv.grad = "ignore"), 
                              data = polity.diffusion.unimputed)
  
  
  unimputed.glm.democ <- glm(democ ~ prop_democ_mimic + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                             family = "binomial", 
                             data = polity.diffusion.unimputed)
  
  unimputed.lm.democ <- glm(democ ~ prop_democ_mimic + democratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                            data = polity.diffusion.unimputed)
  
  
  ## Institutional diffusion of autocracy
  unimputed.re.autoc <- glmer(autoc ~ prop_autoc_mimic + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                              family = binomial(link = "logit"), 
                              control = glmerControl(optimizer = "bobyqa",
                                                     optCtrl = list(maxfun = 5e10))
                              , 
                              data = polity.diffusion.unimputed)
  
  
  unimputed.glm.autoc <- glm(autoc ~ prop_autoc_mimic + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                             family = "binomial", 
                             data = polity.diffusion.unimputed)
  
  unimputed.lm.autoc <- glm(autoc ~ prop_autoc_mimic + autocratic_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                            data = polity.diffusion.unimputed)
  
  ## Institutional diffusion of political participation restrictions
  unimputed.re.rest.part <- glmer(restricted_political_participation ~ prop_rest_pol_part + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                                  family = binomial(link = "logit"), 
                                  control = glmerControl(optimizer = "bobyqa",
                                                         optCtrl = list(maxfun = 5e10))
                                  ,  
                                  data = polity.diffusion.unimputed)
  
  
  unimputed.glm.rest.part <- glm(restricted_political_participation ~ prop_rest_pol_part + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                                 family = "binomial",  
                                 data = polity.diffusion.unimputed)
  
  unimputed.lm.rest.part <- glm(restricted_political_participation ~ prop_rest_pol_part + rest_pol_part_neighbor + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                                data = polity.diffusion.unimputed)
  
  ## Institutional diffusion of free media
  unimputed.re.free.media <- glmer(free_media ~ prop_gov_censorship_mimic + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2) + (1|year) + (1|names), 
                                   family = binomial(link = "logit"), 
                                   control = glmerControl(optimizer = "bobyqa",
                                                          optCtrl = list(maxfun = 5e10))
                                   ,  
                                   data = polity.diffusion.unimputed)
  
  unimputed.glm.free.media <- glm(free_media ~ prop_gov_censorship_mimic + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2),
                                  family = "binomial",
                                  data = polity.diffusion.unimputed)
  
  unimputed.lm.free.media <- glm(free_media ~ prop_gov_censorship_mimic + gov_censorship_neighbor + gle_cgdpc + wdi_gdpcapgr + civil_war  + poly(year_num, 2),
                                 data = polity.diffusion.unimputed)
  
  ## Institutional diffusion of executive constraints
  unimputed.re.exec.const <- glmer(executive_constrainted ~ prop_exec_constraint + exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2) + (1|year) + (1|names), 
                                   family = binomial(link = "logit"), 
                                   control = glmerControl(
                                     optimizer = "bobyqa", optCtrl = list(maxfun = 5e10)
                                   ), 
                                   data = polity.diffusion.unimputed)
  
  unimputed.glm.exec.const <- glm(executive_constrainted ~ prop_exec_constraint + exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                                  family = "binomial", 
                                  data = polity.diffusion.unimputed)
  
  unimputed.lm.exec.const <- glm(executive_constrainted ~ prop_exec_constraint + exec_constraint_neighbor  + gle_cgdpc + wdi_gdpcapgr +  civil_war + poly(year_num, 2), 
                                 data = polity.diffusion.unimputed)
  
})

summary(unimputed.re.democ)
summary(unimputed.re.autoc)
summary(unimputed.re.rest.part)
summary(unimputed.re.free.media)
summary(unimputed.re.exec.const)

summary(unimputed.glm.democ)
summary(unimputed.glm.autoc)
summary(unimputed.glm.rest.part)
summary(unimputed.glm.free.media)
summary(unimputed.glm.exec.const)

summary(unimputed.lm.democ)
summary(unimputed.lm.autoc)
summary(unimputed.lm.rest.part)
summary(unimputed.lm.free.media)
summary(unimputed.lm.exec.const)

save(unimputed.re.democ, file = "model_output/unimputed.re.democ.rda")
save(unimputed.re.autoc, file = "model_output/unimputed.re.autoc.rda")
save(unimputed.re.rest.part, file = "model_output/unimputed.re.rest.part.rda")
save(unimputed.re.free.media, file = "model_output/unimputed.re.free.media.rda")
save(unimputed.re.exec.const, file = "model_output/unimputed.re.exec.const.rda")

save(unimputed.glm.democ, file = "model_output/unimputed.glm.democ.rda")
save(unimputed.glm.autoc, file = "model_output/unimputed.glm.autoc.rda")
save(unimputed.glm.rest.part, file = "model_output/unimputed.glm.rest.part.rda")
save(unimputed.glm.free.media, file = "model_output/unimputed.glm.free.media.rda")
save(unimputed.glm.exec.const, file = "model_output/unimputed.glm.exec.const.rda")

save(unimputed.lm.democ, file = "model_output/unimputed.lm.democ.rda")
save(unimputed.lm.autoc, file = "model_output/unimputed.lm.autoc.rda")
save(unimputed.lm.rest.part, file = "model_output/unimputed.lm.rest.part.rda")
save(unimputed.lm.free.media, file = "model_output/unimputed.lm.free.media.rda")
save(unimputed.lm.exec.const, file = "model_output/unimputed.lm.exec.const.rda")
